library(clusterProfiler)
library(ggplot2)
library(enrichplot)
library(readxl)
library(ggpubr)
library(gridExtra)
library(pheatmap)
library(DESeq2)
library(RColorBrewer)


color = colorRampPalette(rev(brewer.pal(n = 7, name =
                                          "RdYlBu")))(100)

path_barplot = function(all_go, top = 1000000) {
  path_ratio = do.call(rbind, sapply(all_go$GeneRatio, strsplit, '/'))
  path_ratio = apply(path_ratio, 2, as.numeric)
  
  bk_ratio = do.call(rbind, sapply(all_go$BgRatio, strsplit, '/'))
  bk_ratio = apply(bk_ratio, 2, as.numeric)
  
  plot_height = path_ratio[, 1] * bk_ratio[, 2] / path_ratio[, 2] / bk_ratio[, 1]
  
  cols = -log10(all_go$pvalue)
  
  df = data.frame(pathways = all_go$Description,
                  fc = plot_height,
                  'logp' = cols)
  df = df[order(df$logp, decreasing = T), ]
  if (top < length(cols))
    df = df[1:top, ]
  
  df = df[order(df$logp, decreasing = F), ]
  df$pathways = factor(df$pathways, levels = df$pathways)
  temp = df$logp
  temp1 = quantile(df$logp)
  temp[temp > temp1[3] * 5] = max(temp[temp <= temp1[3] * 5])
  temp[temp < temp1[3] / 5] = min(temp[temp >= temp1[3] / 5])
  df$logp = temp
  name_length = max(apply(df, 2, nchar)[, 1])
  ggplot(data = df, aes(x = pathways, y = fc,
                        fill = logp)) +
    geom_bar(stat = 'identity',
             width = 0.9,
             color = 'black') + coord_flip() +
    scale_fill_distiller(palette = "RdYlBu") +
    theme(
      text = element_text(size = 26),
      axis.text = element_text(size = 15 * 50 / name_length),
      legend.text = element_text(size = 15) ,
      panel.background = element_blank(),
      axis.line = element_line(colour = "black"),
    )
}


path_dotplot = function(all_go, top = 1000000) {
  if ('GeneRatio' %in% names(data.frame(all_go))) {
    path_ratio = do.call(rbind, sapply(all_go$GeneRatio, strsplit, '/'))
    path_ratio = apply(path_ratio, 2, as.numeric)
    path_ratio = path_ratio[, 1] / path_ratio[, 2]
    
    
    df = data.frame(
      pathways = all_go$Description,
      GeneRatio = path_ratio,
      Count = all_go$Count,
      logp = -log10(all_go$pvalue)
    )
  }
  
  else {
    path_count = sapply(strsplit(all_go$core_enrichment, '/'), length)
    path_ratio = path_count / all_go$setSize
    df = data.frame(
      pathways = all_go$Description,
      GeneRatio = path_ratio,
      Count = path_count,
      logp = -log10(all_go$p.adjust)
    )
  }
  
  
  
  df = df[order(df$GeneRatio, decreasing = T), ]
  if (top < length(path_ratio))
    df = df[1:top, ]
  
  df = df[order(df$GeneRatio, decreasing = F), ]
  df$pathways = factor(df$pathways, levels = df$pathways)
  
  
  temp = df$logp
  temp1 = quantile(df$logp)
  temp[temp > temp1[3] * 5] = max(temp[temp <= temp1[3] * 5])
  temp[temp < temp1[3] / 5] = min(temp[temp >= temp1[3] / 5])
  df$logp = temp
  name_length = max(apply(df, 2, nchar)[, 1])
  ggplot(data = df, aes(x = pathways, y = GeneRatio)) +
    geom_point(aes(size = Count, color = logp)) + coord_flip() +
    scale_color_distiller(palette = "RdYlBu") + theme(
      text = element_text(size = 26),
      axis.text = element_text(size = 15 * 50 / name_length),
      legend.text = element_text(size =
                                   15)
    )
  
  
}


go_analysis = function(data,name) {

  genes = data
  genome = 'org.Hs.eg.db'
  eg <-
    bitr(genes,
         fromType = "SYMBOL",
         toType = "ENTREZID",
         OrgDb = genome)
gene <- eg$ENTREZID
  
all_go <-
    enrichGO(
      gene = gene,
      OrgDb = genome,
      ont = "BP",
      pAdjustMethod = "BH",
      pvalueCutoff = 1.1,
      qvalueCutoff = 1.1,
      readable = TRUE
    )
  all_go=all_go[all_go$pvalue<0.05,]
  
  
  if (!is.null(all_go)) {
    combined_go <- as.data.frame(all_go)
    combined_go <- combined_go[order(combined_go$p.adj), ]

    if (nrow(all_go) > 1) {
      curheight = 900
      if (nrow(all_go) < 20)
        curheight = 150 + 45 * nrow(all_go)
      
      barname = paste('./figures/sFigure_1BE_',name,
                      '_GO_barplot.pdf',
                      collapse = '',
                      sep = '')
      
      path_barplot(all_go, top = 20)
      ggsave(barname, width = 10, height = curheight / 100)
    
    }
  }
  
}


kegg_analysis = function(data, name) {


  
 genes = data
  genome = 'org.Hs.eg.db'
  eg <-
    bitr(genes,
         fromType = "SYMBOL",
         toType = "ENTREZID",
         OrgDb = genome)
gene <- eg$ENTREZID
  
  kk <-
    enrichKEGG(
      gene,
      organism = 'hsa',
      pvalueCutoff = 1.1,
      pAdjustMethod = "BH",
      qvalueCutoff = 1.1
    )  
  kk=kk[kk$pvalue<0.05,]
  if (!is.null(kk)) {

    if (nrow(kk) > 1) {
      curheight = 900
      if (nrow(kk) < 20)
        curheight = 150 + 45 * nrow(kk)
      barname = paste('./figures/sFigure_1BE_',
                      name,
                      '_KEGG_barplot.pdf',
                      collapse = '',
                      sep = '')
      path_barplot(kk, top = 20)
      ggsave(barname, width = 10, height = curheight / 100)
      

    }
  }
  
}



